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ABSTRACT 


When a radar pulse impinges upon a target, the resultant 
scattering process can be solved as a linear time-invariant 
(LTI) system problem. The system has a transfer function with 
poles and zeros. Previous work has shown that the poles are 
independent of the exciting waveform and target's aspect, but 
they are dependent on the target's structure and geometry. 
This thesis evaluates the resonance estimation performance of 
two signal processing techniques: the Kumaresan-Tufts 
algorithm and the Cadzow-Solomon algorithn. Improvements are 
made to the Cadzow-Solomon algorithn. Both algorithms are 
programmed using MATLAB. Test data used to evaluate these 
algorithms includes synthetic and integral equation generated 
Signals, with and without additive noise, in addition to new 
experimental scattering data from a thin wire, aluminum 


spheres and scale model aircraft. 
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I. INTRODUCTION 


When a radar pulse impinges upon a target, the resultant 
scattering process can be solved as a linear time-invariant 
(LTI) system problem. The system has a transfer function with 
poles and zeros. This description is provided by the 
Singularity Expansion Method (SEM) developed by Dr. Carl Baum 
(Ref. 1], of the Air Force Weapons Laboratory at Kirtland Air 
Force Base. Baum's SEM describes the induced current and 
scattered field using damped sinusoidal waveforms which 
resonate with complex natural frequencies unique to the 
object. These frequencies are determined by the object's 
composition and structural geometry. Natural frequencies are 
also the complex poles of the transfer function. These poles 
are independent of the incident electromagnetic excitation, 
including aspect and polarization, as initially postulated by 
Moffatt and Mains [Ref. 2]. Morgan [Ref. 3] has shown that a 
target scattering response (following illumination) can be 
represented as a weighted expansion of complex natural 
resonances whose poles are independent of the incident 
electromagnetic excitation. The problem of classifying a 
radar target can be solved, in principle, by determining the 


pole positions of the target's response. 


Although this idea is not new, early attempts to 
demonstrate the practicability of such an identification 
system have produced poor results due to the presence of noise 
in the system. Two separate signal processing algorithms have 
been developed by kKumaresan-Tufts [Ref. 4] and Cadzow- 
Solomon [Ref. 5} to locate target poles with a high 
degree of accuracy, in the presence of noise. The Kumaresan- 
Tufts algorithm is intended for purely auto-regressive (AR) 
Signals, while the Cadzow-Solomon algorithm is intended for 
auto-regresSive moving-average (ARMA) signals. This thesis 
evaluates the viability of these two signal processing 
techniques by using new experimental electromagnetic 
scattering data. The use of this method improves 
identification of natural resonances in noisy signals with the 


correct selection of problem parameters. 


A. THE PROBLEM 

The approach employed to classify radar targets is 
comprised of two steps, based on natural resonances. The 
first step locates the position of the poles of the targets of 
interest. Numerical analysis of integral equation techniques 
will derive the poles for simple targets, e.g., a thin wire. 
For more complicated targets such as aircraft, the poles must 
be extracted from actual measurements of the target's response 


to incident electromagnetic excitation. The information 


collected can then be used to form a database for comparison 
to the data obtained in actual field use. 

The second step for target classification is’ the 
comparison of the data obtained in field measurements with the 
data contained in the database. The target is classified 
based on the closest data match. One method of accomplishing 
the classification is to use the same signal processing 
technique that was used to form the database initially. Speed 
is an important factor in the classification process, as the 
time required to achieve classification must be less than the 
time for the target to become a threat. Another, more 
efficient way to perform the comparison is to employ the use 
of annihilation filters, as proposed by Morgan and Dunavin 
(Ref. 6]. This approach employs energy cancellation 
based upon the location of the target's poles. For example, 
a total target response (including early-time and late-time) 
may be subjected to a bank of annihilation filters. The 
filter that corresponds to the proper target will exhibit, in 
theory, a response coincident with the driven portion of the 
target response while annihilating the signal during the late- 
time portion of the response. 

A system used for classification of radar targets would 
require a bank of annihilation filters corresponding to each 
Class of target. The system could then determine the class of 
target, based on the filter producing the minimal amount of 


output energy. 


B. BACKGROUND 

Transient electromagnetic scattering can be described by 
showing an incident field in free space illuminating a 
perfectly conductive finite-sized object. Figure 1 depicts 
induced current on a scattering body. The magnetic field 


integral equation (MFIE) 


UNIT NORMAL 


E1G 4) 


— 
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Figure 1 Transient electromagnetic scattering 


describes induced current on the surface of a scattering body 


as: 


, at = 7 
T(E, t) <2HxHi (F,t)+ff R(z,F/, ) TU, ¢-4Z-21 as (2) 
where, 


¢ J is the electric current density, 


¢ 7 is the unit normal vector outward to the surface, 


- H’ is the incident magnetic field intensity at the 
surface, 


° K is the dyadic Green's function kernel 
and, the principal-value (PV) type integral excludes the point 
r=r’. The cross product, 2nxH*’, represents the physical 
optics portion of the induced current. The principal-value 


surface integral, S,,, represents the contribution due to 


induced physical optics currents other than the contribution 
from the cross product. This current represents "feedback" 


effects from currents at other points on the object. 
When H7=0, Equation (1) solutions become the source-free 


("natural") current modes of the scattering problem. These 


modes may be described by the sum of the exponentially damped 


sinusoids, or: Jd,(r)exp(s,t), where s, are the poles or 


natural resonance frequencies in complex conjugate pairs of 


the form 
S,=O,+j'O,.- (2) 
where, 
° 6, = the damping rate in Nepers/sec, 


° W@W, = the frequency in radians/sec. 


These resonance frequencies are functions of the geometry and 
composition of the scattering object. Although the number of 


poles is, in theory, infinite for any given object, only a 


finite subset will be excited, due to the finite bandwidth of 
the incident field. 
Figure 2 represents a plane wave impulse illumination,in 


which part of the object has been illuminated from the 


IMPULSIVE 
INITIAL 


IMPACT 
POINT 
AT t=@ 


INCIDENT FIELD 
WAVEFRONT 





Figure 2 Plane wave impulse illumination 


incident field wavefront. The scattered field is composed of 
two parts: the early-time driven response and the late-time 
natural mode response. The early-time driven response occurs 
as the excitation wavefront travels over the object and each 
point becomes excited. The late-time natural mode response 
occurs, by definition, when the excitation has been completed 
over the complete body. Throughout the early-time phase, the 
impulsive plane wave incident field will be identically zero 


at all points on the surface except on the conformal ring 


where the intersection with the wavefront occurs. This area 
is indicated by the dotted line in Figure 2. This conformal 
ring changes cross-sectional shape and position as_ the 
wavefront moves. The induced surface current on the ring is 
therefore composed of the physical optics term plus the 
feedback current, as described in Equation (1). Thus, no 
induced current is present at the area of the object ahead of 
the wavefront. 

The back-scattered far-field, resulting from the surface 
current on the object, will be of the form 


ak 


Fs (- = _o 
BN iS Ie 4ncr ot 





= s/ 
[ [exo ez 1) ‘ds’. (3) 
S 


The unit vector D indicates the direction of the plane wave's 
propagation. The back-scattered far-field for a fixed point 
is then a result of integrating the current at every point on 
the object's surface. Thus, the back-scattered far-field will 


be of the form 


H*(-rp, t) =10 eek =n, t) +e ie (-rp, e) exp (s,t)] - (4) 


N=s=-* 
neo 


Two individual terms emerge from these calculations. The 


first term in Equation (4) represents the physical optics 
scattered field generated by the 2nxH’ driven current. The 


second term represents the scattered field produced by the 


source-free wake current or "feedback" current. The physical 
optics scattered field is highly aspect dependent. 

During the early-time portion of the target's response, 
the scattered field produced by the "feedback" current term in 


Equation (4) contains exponential resonance terms with time- 


varying coefficients H,, as described by the Singularity 


Expansion Method (SEM) [Ref. 1]. This form of the SEM 
expansion is termed "Class 2" and is presented analytically by 
Morgan [Ref. 7]. For the monostatic scattering case, 
the transition from early-time to late-time in the target's 


response will occur at 


At=t+(1+cosa)-=, (5) 


after the leading edge of the scattered pulse returns to the 


radar antenna, where, 

* t is the incident pulse width, 

¢ Dis the target's largest dimension, and 

* @ is the aspect of the target. 
The term a also represents the angle between the target's 
largest dimension and the direction of wave propagation, where 
c is the velocity of propagation or speed of light. At this 
transition time, the physical optics field vanishes. During 
the late-time portion of the target's response only the second 


term in Equation (4) remains. However, it is now comprised of 


constant coefficients, H,. This form of the late-time SEM 


expansion is termed "Class 1". 

The early-time scattered field is therefore constituted of 
both a physical optics term and a "Class 2" SEM expansion with 
time-varying coefficients, while the late-time scattered field 
is described by a simple "Class 1" SEM expansion with constant 


coefficients. 


II. POLE EXTRACTION ALGORITHMS 


Obtaining the target response is the first step to be 
completed in a Non-Cooperative Target Recognition (NCTR) 
system. Once the response is obtained in time domain, the 
target's poles may be located. Numerically, these poles may 
be described by the damping rate o and the radian frequency 
® or, in other words, the complex number indicating the pole's 
location in the s-plane or z-plane. The location of the poles 
may also be described analytically by solving the boundary- 
value electromagnetic problem. This may be easily done only 
for simple shape targets like thin wires and spheres. 
However, for general targets, detailed knowledge of the 
target's composition (dimensions and materials) and, also, 
access to a large amount of computing power is required. 

The next section of this thesis discusses the analytical 
method which is used to compare the poles extracted from 
measurements taken with a thin wire. The viability of the 


Cadzow-Solomon algorithm will also be tested. 


A. EARLY METHODS 
The NCTR system is used to discriminate and to classify 


targets with relatively similar shapes and dimensions. In 


10 


such cases, the location of the poles from these targets will 
be relatively close, requiring a high degree of precision in 
measurement and estimation. A basic method in signal 
processing is to obtain the signal's spectrum using the Fast 
Fourier Transform (FFT). The FFT algorithm yields results in 
a short amount of time. However, the FFT may be used only as 
a tool to add information to the basic methodology when 
locating radar target natural resonances. This is because the 
frequency resolution, equal to the reciprocal of the time 
interval between the samples, is of the order of several MHz 
or higher. Another limitation of the FFT is that it provides 
only real frequencies within a signal, while both the damping 
rate and the frequency of the poles are required for target 
classification. 
Thus, other methods providing more accurate results needed 
to be developed. 
1. Direct Minimization 
As previously described, the transient scattered 


Signal waveform will have the form 


y(t) =y,(t) (u(t) -u(t-7,)] ty, (t) ‘u(t-T,) +N(t) (6) 


where, 


° yz(t) is the early-time portion, 
° y,(t) is the late-time portion and 


« N(t) is the undesirable noise and clutter. 
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Morgan [Ref. 7] presents the way to determine the natural 
resonance parameters by modeling the late-time scattering 


response as a sum of damped sinusoids, 


¥,(t) =)> A,exp (0, t) cos (w,t+,) (7) 


n=1 


where y represents the modeled signal and 4A,‘exp(j®@,) 


represents the aspect dependent residues. The digital domain 


form of the model is 


N 
¥(k:At) =y,=)> A,exp (0,°kA t) cos (w,"kA t+9,) (8) 


n=1 


¢ A, 1s the amplitude, 
* 6, is the damping rate, 
* @, 1s the radial frequency, and 
¢ @, is the phase. 
After comparing the modeled signal y, to the actual received 


discrete signal y,, the mean square error may be obtained as 
ef=(y,-y,)?. The four parameters of the model (A,, $,, 9,, ®,) 


may then be adjusted in order to derive the mean square error 
minimum. As this is a multi-dimensional, highly non-linear 


Minimization problem, it is computationally inefficient. 
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2. Prony's Method 
Blaricum and Mittra [Ref. 8] present a novel 
approach for systematically deriving the complex poles and 
residues of a target from a set of time domain data. The 
method is based on Prony's algorithm, used to model the late- 
time response of a radar target. The set of time domain data 
is the discrete set of sampled transient values of the impulse 


response I(t,) or I,.- The method is based on the fact that 


I, must satisfy a difference equation of order N, written as 


N-1 
S Qp tout Iek ’ (9) 
p=0 

where p+k=0,1,...,2N-1 and 2N are the data samples used. This 


leads to the solution of a matrix equation 


ee I 
- . " ay Ie 
ees Zea} [yal [Inez 
: (10) 
ay lon 
Ly tyes Ioy-1 


After the coefficients a, are found, the method solves the 


roots of a polynomial as 


N N- = 
z *-a,z"1-..-a,=0. (11) 


If d, (m=1,2,...,N) represents those roots, then the poles of 
the system model in the z-plane are the roots d,. The poles 


in the s-plane may be found using the formula 
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(12) 


where At is the time-stepping interval used in obtaining the 
sampled data. The matrix in Equation (10) is in the form of 
a transposed Vandermonde matrix, whose inverse can be computed 
in closed forn. This method requires at least 2N equally 
Spaced transient data samples to find N poles. If greater 
than 2N samples are desired, then one may obtain a least 
squares type fit to the matrix in Equation (10). Since the 
system order is not known a priori in the NCTR problen, 
Blaricum and Mittra [Ref. 9] present two schemes’ for 
determining the number of poles contained in the transient 
response. The first is the Householder orthogonalization 
method, and the second is an eigenvalue method. When the data 
are noisy, Blaricum and Mittra [Ref. 9] indicate that the 
eigenvalue method is the better method. This method will be 
discussed in Section 3 of this chapter when a bias 
compensation method is examined. 
3. Kumaresan-Tufts Algorithm 

Kumaresan and Tufts (Ref. 4] modified Prony's method 
by using the "backward linear prediction" technique and 
"Singular value decomposition" (SVD) to alleviate the 
sensitivity to noise and the need for a priori knowledge of 


the system order. 
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Analytically, Kumaresan-Tufts algorithm modifies 
Prony's method as follows: 
¢« The system order is intentionally overestimated. The 
excess order provides the flexibility to the system to 
model the noise, improving the estimation of parameters of 
exponentially damped sinusoidal signals in noise. 


¢ Singular value decomposition alleviates severe ill- 
conditioning of the data matrix. 


* The causality of the system is used to separate the 
computed signal's poles from the spurious computed noise 
poles introduced by the overestimated system order. 

The separation of the signal poles from the noise 
poles is the result of the "backward prediction" that causes 
the signal poles to fall outside the unit circle in the z- 
plane, while the extraneous "noise poles" fall inside the unit 
Gircle. 

Moderately large values of system order are essential 
in improving the accuracy of the pole location estimates. 

a. System Model 

Suppose that the N samples of the observed time 
domain data of a response signal y(n) consists of samples of 


M exponentially damped signals in white Gaussian noise w(n), 


as represented by 


M 
y(n) => ayexp(s,n) +w(n), n=0,1,...,N-1. (13) 
k=1 


The following linear prediction equations can be formed, using 


real time domain data in the backward direction. 
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Yley 0 YUST VA EE a) y(t) 
y(3) yl4) .. y(L+2) a(2)}__] y(2) (14) 


i a a(L) (N-L) 
(N-L+1) .. .. y(N) 
The terms of the above equation may be also represented by 
D:a=-h (14a) 
where, 
- Dis the data matrix (N-L)xL, 


°° ais the vector of backward prediction coefficients 
(LxX1), and 


¢ h the data vector (N-L)x1. 
In the above matrix equation, L represents the overestimated 
system order, chosen to satisfy the inequality MSLSN-M. The 
matrix equation is solved with the SVD method, using the 
pseudoinverse matrix D*, as a=D*:(-h). Here the coefficients 


a, correspond to those in Equation (11), where the roots of 


the polynomial define the poles of the system model in the z- 
plane. 

As with Prony's method, the Kumaresan-Tufts method 
is intended for purely auto-regressive (AR) signals. Both 
methods, therefore, use the late-time portion of the response 


Signal when addressing the target classification problem. 
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b. Singular Value Decomposition 


Singular value decomposition (SVD) is the basic 
technique on which the Kumaresan-Tufts algorithm is based. 
There are two main advantages with SVD: 


¢ It helps to alleviate the effects of ill-conditioning of 
the data matrix. 


¢ It separates the signal poles from the noise poles. 
The following discussion of SVD is a synopsis of material 
taken from Strang [Ref. 10]. 


The SVD is closely associated with the eigenvalue- 
eigenvector factorization of a symmetric matrix : A=QA‘Q7. 


SVD makes a Similar factorization, but for any (mxn) matrix 


A, as follows: 


A-One. (15) 


where Q,, Q, are two orthogonal matrices and 2 is the diagonal 


(but rectangular) matrix with its positive entries (also 


called sigma, in the form of 0o,,0,,...,6,). These entries are 


the singular values of A, filling the first r places on the 
main diagonal of 2, where r is the rank of A. These entries 
Sigma are also the square roots of the nonzero eigenvalues of 


both AA? and A’A. The columns of Q, (mxm) are eigenvectors of 


AAT, and the columns of Q, (nxn) are eigenvectors of A’7A. 
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The SVD works well for numerically stable 


computations. An initial reason is that Q, and Q, are 


orthogonal matrices that never change the length of a vector. 
Since ||Qx||?=x 70 70x=||x||?, multiplication by Q cannot destroy 
the scaling. A second reason is that SVD gives a more stable 
measure of the rank of A. 

The prime use of SVD is to solve every linear 


system in the form of Ax=b. For every matrix A (m xX n), which 


can be factored into A=0,‘D:0,’, another matrix can be defined. 


This matrix will be the pseudoinverse of A, as follows: 


A= Oo Ole (16) 


where, Q,, 0, are the same orthogonal matrices found with the 


SVD. The value %* is the diagonal (n x m) matrix, with its 





entries on the main diagonal being — a nr a The 
0, 9, 0, 
optimal solution of Ax=b is 
x*=A*b (17) 


that is, the minimum length solution to the nearest equation, 


Ax=p. As the column space of A* is in the row space of A, 


then x* is always in the row space of A. 
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To solve the system of equations in Equation (14), 
the pseudoinverse of matrix D may be found as in Equation (16) 
in the form D*=V~*:UT, where Z* will be a (L xX (N-L)) matrix, 
whose L singular values are the reciprocals of those found in 
the & matrix. The minimum length least squares solution of 
Dia=-h is a=D**(-h) =VE*:U"(-h). 

c. Bias Compensation 

In the discussion of the SVD above, 2 is a (N-L)xXL 
Giagonal matrix, with its entries being the square roots of 
the nonzero eigenvalues of both DD? and D’D. In the noiseless 


case, the prediction equations are satisfied exactly. If the 


system order has been overestimated as L, when M is the actual 
order of the system, the diagonal of XZ splits into a signal 


subspace with M positive singular values of 6o,,0,,...,06, and 


a subspace with L-M zero values. For the noisy signal case, 
the first M positive signal values are perturbed into a noisy 


Signal subspace 
O,,05,+++,0, With eigenvalues \j2zAj2...2Ay, (18) 


and a noise subspace with L-M values 


e 6 / / / 
Omer Ouse +e, O0, With eigenvalues Amy2Ayi22...2A,. 019) 


Kumaresan and Tufts noticed [Ref. 4] that the noise 


perturbed the signal's singular values, the reason for the 
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bias towards the unit circle in the z-plane in the pole 
position estimates. Kumaresan and Tufts proposed a 
compensation method which moves the poles back towards the 
center of the z-plane. This method is based on averaging the 


smallest L-M singular values, Oy,;,Om0----,90,, and then 


subtracting this average value from each of the first M 


Singular values, 0,,0,,...,0,- The smallest L-M values are 


then set to zero and a new matrix, 2, is used to compute the 
pseudoinverse D*. Although Kumaresan and Tufts claim that 
this method gives better results, the testing completed for 
this thesis shows that the bias towards the center of the z- 
plane is greater than the bias error that the noise makes. 
This result causes much more perturbation on the pole position 
estimates. 

Norton fRef. 11] proposed another 
compensation method based on the Eigenvalue Shifting Theoren. 
The noisy data matrix may be described by D=S+N, where S is 
the noiseless signal data matrix and N is the noise data 
matrix of the wide-sense stationary white noise process, v,, 


written as 


Vay, 


(20) 


20 


The expected value of DD’ may be found from 


DD T=E[(S+N) (S+N)7] =E(SS 7) +E[SN™+E[NS7™]+E[NN™]. (21) 


As white noise has a mean of zero, the two terms E[SN7] and 
E(NS*] become zero. Since the noise is wide-sense stationary, 


E(NN7]=0.:I, where o, is the noise variance and I is the 


identity matrix. As S is deterministic, E[SS7)=SS7. This 


leads to the formula 


E[DD7] =SS T+02-I. (22) 


The eigenvalue shifting theorem [{Ref. 11] describes the case 


when the eigenvalues of SS*7 are 4,;, where the eigenvalues of 


E[DD"] are A, +0. Therefore, the eigenvalues of DD‘, which 


are the squares of the singular values of D, are increased by 


the noise variance, o°. 


In the noiseless case, these singular values would 


be zero. If the noise singular values are squared and 


averaged, an unbiased estimation of of may be obtained. 


Norton (Ref. 11] proposed correcting the signal singular 


values by first subtracting the noise variance o% from the 


squares of the first M diagonal entries of matrix 2 and then 


taking the square root of the reduced values as the new 
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diagonal entries of matrix Z (while the noise singular values 


are set to zero). As in the Kumaresan-Tufts bias compensation 
method, the pseudoinverse D* can be found from the new matrix 


2ive 

Although Norton's method seems to be a more correct 
bias compensation method, both methods are based on the fact 
that a decision has to be made about the actual order of the 
system. The separation between the signal eigenvalues and the 


noise eigenvalues is not readily obvious, because the 
eigenvalues of matrix DD’? or D"D appear to steadily decrease. 


ad. Performance and Earlier Results 

Kumaresan and Tufts tested the algorithm, obtaining 
good results using synthetic data with various levels of white 
Gaussian noise, as low as 15 dB SNR. Norton [{Ref. 11] 
developed the algorithm as a computer subroutine and tested it 
with various sets and types of data. Synthetic data were 
tested at various SNR values, ranging from 90 dB to 7 GB. 
When the SNR decreased, the error in pole position estimates 
increased. Norton claimed that the algorithm gave good 
results for SNR220 dB. Norton also used simulations of thin 
wires produced by Morgan's Time-Domain thin wire integral 
equation (TDIE) computer program to test the algorithm. The 
results of those tests illustrated the aspect independence on 


estimated poles of the target, [Ref. 11]. 
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Larison (Ref. 12] programmed the - algorithm 
using Fortran and tested synthetic and thin wire integral 
equation data at various SNR's ranging from 90 dB to 7 GB. 
Like Norton, he claimed good results in extracting the low 
frequency poles position with SNR as low as 20 dB by using 
Norton's bias compensation method. lLarison maintained that 
the two most critical parameters are: 


¢ to select the appropriate starting point to begin the data 
processing 


¢ to select the appropriate system order to provide the best 
possible results. 


The algorithm has been tested by the author of this 
thesis, after developing the algorithm using MATLAB. The 
algorithm was tested using synthetic data at various SNR's 
ranging from 30 adaB to 10 dB, without using the bias 
compensation method. The synthetic signal response is based 
on ten pole pairs within a frequency range of 1-10 GHz, with 
a medium Q damping factor using k=0.7. This chosen value of 
k simulates typical expected levels of damping from measured 
scattering responses of real targets (e.g., the thin wire). 
Damping factor and SNR level are discussed, in detail, in 
Chapter III, Section A, of this thesis. Table I shows the 
poles and residues used for the testing. 

Figures 3 to 10 illustrate the evaluation of the 
Kumaresan-Tufts algorithm. Figures 3 to 6 illustrate the 
synthetic signal waveform for various SNR's. Figures 7 and 8 


illustrate the spectrum of the synthetic signal for SNR=30 dB 
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and SNR=10 aB. Figures 9 and 10 illustrate the true and the 
algorithm poles in both the s-plane and z-plane, for various 
SNR's. For the above simulation, the values used for N and L 
were 300 and 60, respectively, while M (actual order) was 20. 
The simulation revealed that as the noise increases relative 
to the signal, there is a bias of all the poles towards the 


unit circle, specially the higher frequency poles. 


Table I MEDIUM Q SYNTHETIC POLES 


aS 


=e ea a 
=o ale ae 
ar [eee 
we See 
Bea to 
—s 
= 
eer 
i 
3 





a 





[GHz] 





r — = = —— 





=20s 73 37.6511 
-2.4935 43.9263 


-2.8498 50.2015 


-~3.2060 56.4767 





=3.5622 62.7518 


24 





amplitude 


amplitude 


0.8 


0.6 


0.4 


0.2 


0 50 100 150 200 250 300 350 400 450 500 


time points 
Figure 3 Generated synthetic signal without noise 
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Figure 4 Generated synthetic signal with 30 dB SNR 
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Figure 6 
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Generated synthetic signal with 20 dB SNR 
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Generated synthetic signal with 10 dB SNR 


26 





500 





500 


0.025 


Amplitude 


0.03 


0.01 





0 Z 4 6 8 10 12 14 16 18 20 


Frequency in GHz 


Figure 7 Spectrum of the SNR=30 dB synthetic signal 
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Figure 8 Spectrum of the SNR=10 daB synthetic signal 
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Figure 9 Kumaresan-Tufts poles in the s-plane 
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Figure 10 Kumaresan-Tufts poles in the z-plane 


28 


B. CADZOW-SOLOMON ALGORITHM 

Many investigators have used both excitation and response 
data to identify linear systems. Cadzow and Solomon [Ref. 5) 
proposed an identification procedure in which both the 
excitation and response are contaminated by white noise. This 
method is based upon the null space characterization of an 
associated "data matrix". The Cadzow-Solomon algorithm is 
used in this thesis to simulate and demonstrate the 
identification problen. 

1. System Model 

The Cadzow-Solomon algorithm is based upon the Auto- 

Regressive Moving-Average (ARMA) model. In an ideal modeling 
Situation, the assumption may be made that the excitation 
Signal x(n) and the response signal y(n) are perfectly 
related by means of an ARMA relationship, such that 


L K 
y(n) =) ajy(n-i) +)" b,x(n-i), (23) 
2=0 


i=1 


as associated with the transfer function 


at -K 
7 J Op snl OR ame 3 Oe 


igi hs), 
1-a,z-...-a,z™ 


(24) 


where, 


¢ a, are the coefficients which correspond to the poles of 
the transfer function, and 
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- b, are the coefficients which correspond to the zeros of 


the transfer function. 


The order of the numerator of the system's transfer function 


is K and the order of the denominator is L. The classic 


modeling problem is to identify the system's coefficients, a, 


and b,, from a finite set of observations of the excitation, 


x(n), and response, y(n), time series. 


Cadzow-Solomon present the set of equations for the 


algorithm in matrix form as follows: 


y(0O) y(-1) .. yl-p)]74 7 [x(0) x(-1) 
y(1)  y(0) .. y(1-p)} ja,|_jx(1) x0) 


(N) y(N-1) ... y(N-p)| l?P) Ix(N) x(N-1) 


or, more concisely 
Y¥pap=Xqby 


where, 


. X(N-q) 


x(-q@) ][p 


. x(1-9)| |b. } (25) 


b 


(26) 


° X, is the (N+1)x(qt+t1) "excitation data matrix", 


p is the (N+1)x(p+1) "response data matrix", 


* a, is the (p+1)x1l "“auto-regressive parameter vector" and 


° by is the (q+1)x1 "moving-average parameter vector". 


The least squares solution for a causal, real data, 


overestimated ARMA model (p2p,qeq) can be found by taking a 
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linear algebraic approach based on the eigenvalue-eigenvector 


decomposition of the data matrix D, ,=[Y,:!-X,] as 
Dp, g’Dp, q’Up=Ag Uz: 1SKsp+q+2. (27) 


The parameter vector is found from the procedure to be 


p+q+2 Ju, (1) an p+q+2 (1) — (28) 


ef 2; 





kK 
kel Ax k=1 Ax 


where u, is the eigenvector associated with the eigenvalue jA,. 


If the actual system order is (p,q), then the system order is 
less than the overestimated order Cp,.G)* At least 
Ss=l1+min(p-p,q-q) of the eigenvalues in the decomposition, as 


described in Equation (27), must be identically zero for the 
noise-free case. 

Following the technique of "backward prediction" as 
used in the Kumaresan-Tufts algorithm, the matrix equation 


(Equation 25) may be modified as follows: 


y(2). y(Le1)) x(a). Ke) 1 
ved) 
bs),  Y(L+2)) x(2) .. x(K+2) | Jaz) | y(2) (29) 
:. " : 0 : 
(N-L+1) .. y(N) |k(N-L) .. x(N-L+K) ||, OSE) 
Lk) 
or, in matrix notation 
a 
[D,.]°|--l=y. (30) 
aA 
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This matrix is associated with the coefficients a, and b, as 


shown in Equations (23) and (24). The size of the excitation 
and response data matrices are (N-L)x(K+1) and (N-L)XL, 
respectively, while the size of the ARMA parameter vector is 
(L+K+1) x1. 

As in the kKumaresan-Tufts case, singular value 
decomposition can also be used in Cadzow-Solomon's algorithm 
to provide the minimum-norm solution. Its use reduces ill- 


conditioning of the data matrix, D 


yxr and separates the signal 


poles from the noise poles across the unit circle. The 
optimal solution of Equation (30) is 
a 


b 








=[Dyx]"+y. oa 


2. Bias Compensation 


When the system's order, (p,q), exceeds the true order 
of a noise contaminated system, (p,q), for p>p and @q, 


Cadzow and Solomon maintain that there will be 


s=1+min (p-p, q-q) (32) 


eigenvalues, which will asymptotically approach (N+1): 0. for 


large N, where o% is the noise variance. The eigenvalues in 
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the noise-free case are zero, and therefore, the parameter 


vector may be found from the equation 


ap 





s “1s 
> gc) P| ‘yu, (1) 2, (33) 
k=0 


(St) bg £6 


where the u, terms are the eigenvectors associated with the 


smallest (s) eigenvalues, A,. 


However, the above method may not be applied exactly 
as described above because 


* the excitation noise, o,, is zero for the radar target 
classification problem, and 


* the q-q, cannot be directly determined since the 
extraneous zeros cannot be separated from the signal zeros 
in the same manner as the poles are. 

Another problem occuring in the Cadzow-Solomon's 
algorithm is that additive noise is different for input and 


output data. Norton [Ref. 11] described that if the input 


data noise is v, and the output data noise is w,, the noisy 


data matrix may be expressed as 
Dyx=[Dy i D,J=Syx+Nyx (34) 
where, N,,=(N,iN,], or in matrix form 
Ve eV 


i =| (35) 


Vy-1 eee Vaex-L 
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We eee Wri 


y - ~ : (36) 
Wyre eoe Wy 


Therefore, the correlation product will be 
EID, . Dyx|=S,x Syx+ HN, Nyx| (37) 


As the additive noise is different for the input and output 


data, the noise correlation matrix E(N,, Nyx] is not of the form 


Co So, the eigenvalue compensation method that Norton 


described as being applicable in the Kumaresan-Tufts algorithm 
is not applicable in the Cadzow-Solomon algorithn. 


By examining the diagonal entries of the data matrix 


product, Dog Dom it may be seen that the first (g) entries 


are close to 1. In the noise-free case, these first (g) 


Giagonal entries are exactly equal to 1. This is a direct 
consequence of the fact that the diagonal entries of D,,'D,, 


contain bias errors proportional to the noise variances o% and 


hee Since the excitation data used in radar target 


classification problems are noise-free, these bias errors are 


proportional only to o%. The bias compensation may be 


obtained by setting the first (g) diagonal entries equal to 1 


and then finding the eigenvalues and eigenvectors from the 
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corrected matrix, D,.g'D,,,~- The parameter vector may then be 


found from 


=P 


b 





u, (38) 


kel A, 


=i 
peg? Aaa Pee? (1) 
kel A, 





3. Earlier Results 

Norton [{Ref. 11] programmed the Cadzow-Solomon 
algorithm in Pascal and tested it on various types of data. 
Using synthetic data, Norton tested the algorithm using three 
SNR values, 30 GB, 20 @B, and 10 @B. The results were good at 
the 30 dB level, but as the SNR decreased, the error in the 
pole position estimates increased. Norton also used 
simulations of thin wire scattering produced by Morgan's TDIE 
computer program to test the algorithm, obtaining better 
results than those obtained from the Kumaresan-Tufts 
algorithm. These results occured because the Cadzow-Solomon 
algorithm takes into account only the input and output data 
and does not require purely late-time data. 

Larison {Ref. 12] programmed the Cadzow-Solomon 
algorithm in Fortran and tested it using synthetic data and 
TDIE data, at various SNR's, ranging from 90 dB to 7 GB. 
Using bias compensation, Larison claimed good results in 
extracting the low frequency pole positions with SNR as low as 


20 GB. 
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Murphy {Refs 33] tested the Cadzow-Solomon 
algorithm using Larison's Fortran programs, using synthetic 
data as well as thin wire integral equation generated data and 
measured data. Murphy generated three separate signals, each 
based on ten pole pairs covering the frequency range of 1-10 
GHz. Depending on the level of damping for each signal, 
Murphy obtained an average error between the true and the 
extracted poles having values on the order of 10% for the 
values 

¢« SNR=10 dB and for HIGH Q, 

« SNR=20 dB and for MEDIUM Q, and 

¢« SNR=30 dB and for LOW Q. 

Murphy made the observations that, 

¢ The best results were obtained by choosing a starting 
point located within several points of the zero crossing 
nearest to the first obvious response of the excitation. 

¢ The best results were obtained by using a data matrix in 
which the overestimated number of poles to the true number 
of poles was of the order of 2.5, while the number of 


asking zeros was equal to the number of true zeros. 


¢ The best results were not always obtained when using the 
bias compensation scheme as proposed by Norton. 


When processing the thin wire data, Murphy calculated the 
feed-forward order of the system by agmatine the length of 
early-time as equal to 2L/c. Recent research at NPS, that is 
as yet unpublished, indicates that this method is not correct, 
as each excited point along the target eroieee its adjacent 
point. This early-time may be described as (1+cosa):L/c, as 


shown in Chapter III, where a is the aspect of the target. 
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When processing the TDIE data, Murphy did not obtain good 
results for SNR=20 dB, in spite of the fact that the poles 
extracted from thin wire measured responses appeared to give 


good results for the low frequency poles. 
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IIIT. ALGORITHM TESTING 


The initial objective of this thesis was to evaluate the 
Viability of the kKumaresan-Tufts and the Cadzow-Solomon 
algorithms. In a radar target classification problem the 
finite duration excitation signal produces both early-time and 
late-time response signals. In this situation, the Cadzow- 
Solomon algorithm is more significant. Thus, the main effort 
of this thesis was changed to evaluate and improve the Cadzow- 
Solomon algorithm using both synthetic and real data. The 
Kumaresan-Tufts algorithm was evaluated and tested only with 
synthetic data, as presented in Chapter II. 

The Cadzow-Solomon algorithm was tested in two phases. 
The first phase of testing used synthetic data, while the 
second phase was performed with thin wire measurement data. 
The synthetic data testing phase attempted to simulate the 
conditions expected from the response of a simple target 
during the presence of a stationary white noise. The thin 
wire data testing phase attempted to evaluate the conditions 
appearing from a real target response. Those conditions were 
then compared with those simulated from the computed Time 


Domain Integral Equation (TDIE) thin wire response. 
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A. SYNTHETIC SIGNAL MODEL 

As the Cadzow-Solomon algorithm is based on the ARMA 
model, the representation in Equation (23) has been used to 
produce the synthetic signal response. 

This Signal response is based on ten pole pairs within a 
frequency range of 1-10 GHz, with a medium Q damping factor 
using k=0.7. The poles were developed in accordance with the 


following equation: 


wn =_1 ink) (39) 
om 


n 


| 


Table II lists the s-plane poles used in the synthetic signal. 


Table II MEDIUM Q SYNTHETIC POLES 


On 
[Grad/s] 


Qh 
aoe) 
N 


6.2752 
12.5504 
18.8256 
25.1007 
31.3759 
37.6511 
43.9263 
50.2015 
56.4767 
62.7518 


1 
2 
3 
4 
5 
6 
v 
8 
9 
0 


al 





The chosen value of k simulates typical expected levels of 
damping from measured scattering responses of real targets, 
(e.g., the thin wire and scale model aircraft targets). The 
sampling frequency used to convert the s-plane poles to z- 
plane poles was 51 GHz, based on N=256 samples over a time 


window of t,=5 nsec: 
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(40) 


1. Coefficient Generator and Recursive Signal Generator 
The coefficients a, are obtained by multiplying the 


h 


terms (Z-2Z,) (2-2,)...(zZ-Z.), where Zz; is the i z-plane pole 


associated with the s-plane pole in the relationship 


Z,=exp[(o,+jo,)°At), (41) 
and equal to the product of those terms with the polynomial 
20a 2 a eee ee (42) 


The coefficients b; are obtained by the inverse partial 


R; th 





fraction expansion using the term , where R, is the i 


ney 
residue of amplitude value 1 and phase difference O for each 
of the signal poles. The time domain signal response of the 
ARMA model is generated via Equation (23) with the values L=20 
and K=19. This procedure has been developed in the MATLAB 
computer program ARMAGEN1.M. The program Wecing appears in 


Appendix A. 
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2. Double Gaussian Smoothing Function Generator 
The excitation signal chosen for generating the test 


Signal response was the double Gaussian waveform, via the 


equation 
br) YW: 9 00 Car Pal ota Bey =>. 90 at Pu oan jy 
with 
ee eee , ,20.147 nsec (44) 
51 
cycle Se ,  1,=0.314 nsec (45) 
T2 
where, [ze 
A= ’ A. =A,-1 ° (4 6) 
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This waveform is a wide Gaussian pulse with a ten percent 


width of T, nsec subtracted from a narrow Gaussian pulse with 
a ten percent width of T, nsec. This results in a bandwidth 


of 1 to 10 GHz. Figure 12 illustrates the double Gaussian 
waveform and Figure 13 illustrates the spectrum of the 
waveform. This procedure has been developed in the computer 
program EXCGEN.M in MATLAB. The program recing appears in 


Appendix B. 
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3. Synthetic Noise Generator 
Synthetic noise was generated to contaminate the 
signal response by adding a time series noise signal to the 
Signal response. The noise was assumed to be wide sense 
stationary and white. To produce a Gaussian distribution, a 
normal distribution function was multiplied with a standard 


deviation value o, computed via the equation 


N 
Sleek 


k=1 (47) 


o*°= variance = re 
N SNR 


where, 

¢ y(n) is the signal response data, 

- N is the number of time points and 

¢ SNR is the Signal-to-Noise Ratio in dB. 
This procedure was developed in the computer program 
NOISEGEN.M in MATLAB. The program listing appears in 
Appendix C. 

4. Spectrum Estimation 
Estimation of the power spectral density (PSD), also 

Called the spectrum of the sampled | response, is 


obtained employing the Fast Fourier Transform (FFT). This 


spectrum may be computed via the equation 


ak 2 48 
S(k) ay LY CK) | (48) 
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where, 

¢« N is the number of time points (power of 2), 

* Y(k) is the FFT of the signal y(t) and 

¢- S(k) is the periodogram spectral estimate of y(t). 
This procedure was developed through the computer program 
SPECTRUM.M in MATLAB. The program listing appears in 


Appendix D. 


B. SYNTHETIC SIGNAL TESTING RESULTS 

This procedure shows the weaknesses of the Cadzow-Solomon 
algorithm by using no bias compensation. Cases were examined 
for the three different SNR's of 30, 20 and 10 GB. 
Overestimation varied from 2:1 up to 5:1, e.g., for 20 true 
poles and 60 asking poles, the ratio is 3:1. The interval 
processed was composed of 200 points, starting where the 
excitation began. 

Figures 11-23 illustrates the results of this effort. 
Although the generated SNR's were 30, 20 and 10 dB, the 
resulting SNR's over the processed window were 2 dB higher. 
Observations made led to the following factors, thus improving 
results as previously obtained from the Cadzow-Solomon 
algorithm. 

* The most significant factor is to select the excitation 


starting point as the point to start the algorithm 
processing. 
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* The second significant factor is to select the system 
order or number of asKing poles. In the case of 20 true 
poles, 60 asking poles gives the most accurate results for 
all SNR‘'s employed. 
- A third significant factor is to determine the number of 
unknown zeros. For synthetic data, the position of the 
low frequency poles was obtained with better accuracy when 
this number was equal to the number of asking poles than 
with the case where the number of asking zeros was equal 
to the number of true poles. 
Figure 23 illustrates the case of K+1=20 (as compared to 
Figure 19 for K+1=60), where (K+1) represents the number of 
asking zeros. In the case where K+1=20, more poles were 
obtained within the unit circle (as compared to the case where 
K+1=60) at the frequencies of the true poles, but with a 
different damping factor. The experimental data indicates the 
following result: 
- There is a bias of all the poles towards the unit circle 
that results in the loss of some poles. The bias is so 
large that it forces the poles to appear on the other side 
of the unit circle as noise poles. 
By examining the spectra of the synthetic signal at different 
noise levels, as illustrated in Figures 15-16, it may be 
noticed that the high frequency poles cannot be completely 
separated from the white noise spectrum when increasing the 
noise level. This is due to the fact that these high 
frequency poles carry very little energy. 

The set of equations in matrix form developed in Equation 


(29) have been developed in the computer program CADSOL1.M in 


MATLAB. The program listing appears in Appendix E. 
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Figure 12 Spectrum of the Double-Gaussian Function 
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Figure 13 Synthetic signal without noise 
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Figure 14 Synthetic Signal with SNR=10 dB 
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Figure 16 Spectrum of the synthetic signal with SNR=10 dB 
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Poles in z-plane, SNR=32 dB, for synthetic signal 
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Figure 19 Poles in s-plane, SNR=22 dB, for synthetic signal 
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Figure 20 Poles in z-plane, SNR=22 dB, for synthetic signal 
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C. THIN WIRE SIGNAL TESTING 
The Cadzow-Solomon algorithm was also tested using thin 
wire measured scattering data. The results have been compared 
with the results obtained using time domain integral equation 
(TDIE) thin wire data. The poles obtained by processing the 
TDIE computed data were assumed to be correct, even though the 
program that computes the currents on the thin wire does not 
take into account the capacitance at the ends of the wire. 
1. Thin Wire Integral Equation Computed Data 
For the thin wire, the natural resonances may be 
determined by solving the integral equation that describesthe 
current flowing on the wire. The result of this type of 
Simulation is the response of the wire to a specified 
excitation field. Morgan [Ref. 14] developed a time- 
domain thin wire integral equation computer routine, based on 
the formulations of Sayre and Harrington [Ref. 15]. 
The wire used for this simulation had a length of 0.1 meter 
and a radius of 1.18 mm. The back-scattering response was 
computed at three different incident aspects, ranging from 30 
to 90 degrees, with 90 degrees representing a broadside 
aspect. The excitation waveform used was the same double 
Gaussian as that used with the synthetic data, as illustrated 
in Figure 11. Figure 30 illustrates the position pole 
estimates obtained using Cadzow-Solomon algorithm and shows 


the aspect independence of the poles for three cases. The 
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five low frequency poles appeared exactly at the same position 
for all aspects. It should be noted that these results, as 
illustrated in Figure 30, appeared to be exactly the same, 
irrespective of any variations in the parameters (No. of 
points - No. of asking poles - No. of asking zeros) used in 
processing the signal. Note that for broadside illumination, 
only the odd-numbered poles appear. Because of the physical 
symmetry of both the thin wire and incident field, the even- 
numbered modes are prevented from being excited by the 
incident field. With illumination at 45 degrees aspect, a 
spectrum with adequate energy only within the bandwidth from 
1 to 8 GHz was produced, as expected. At higher frequencies, 
backscattering is suppressed because most of the energy is 
reradiated near to the specular scattering directions. 
Therefore, only the five low frequency poles are accurately 
obtained for this case. 

The TDIE program generated a response at 30 degrees 
aspect consisting of 255 points over 5 nsec, resulting ina 
sampling frequency of 50.8 GHz. In the case of 45 and 90 
degrees aspect consisting of 240 points over 5 nsec, a 
sampling frequency of 47.8 GHz resulted. These sampling 
frequencies differ from the sampling frequency of 51 GHz used 
in the measured data. 

When processing the TDIE thin wire eaey points from 
160 to 200 were used to run the algorithm starting at the 


point where the excitation starts. The number of asking poles 
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ranged from 28 (for the 90 degrees aspect) to 40 (for the 30 
degrees aspect). The number of asking zeros or feed-forward 
order of the system was either the same as the number of 
asking poles or, was calculated by determining the length of 


early-time. This early-time was computed from the formula 
ders 

t.=—°(1+cosa) (49) 
C 


where, 

¢ Lis the length of the wire and 

* @ is the aspect angle from end-on orientation. 
This value of time was then converted to the appropriate 
number of time points, as based on the sampling time interval 


Of 


E 
n,= integer|-£2]+1 (50) 


where, 
* T is the sampling interval, and 


* mn, is the feed-forward order of the system. 
This value of n, is the minimum number of asking poles that 


may be used, presenting the number of delays in the z- 
transform for the early-time of the system. 

Another consideration in the processing of these data 
was the scaling of the excitation waveform and its position 
with respect to the computed response. Although the data were 


generated using a double Gaussian waveform with a 1 Volt peak 
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amplitude, the excitation waveform had approximately the same 
peak amplitude as the response waveform. Such scaling does 
not change the frequency contents of the exciting waveform and 
gives better results in the pole extraction due _ to 
minimization of some of the effects of ill-conditioning in the 
data matrix. 

To position the driving waveform with respect to the 
response waveform, the time difference between the excitation 
of the first and last point of the wire may be computed. This 


time interval is represented as 
ag ke 
Cuelay= ~~ 'COSa, (51) 


where a is the aspect angle. The maximum absolute value of 
the response occurs when the last point of the wire is being 
excited. The information derived from the time interval, 


Cuctayy and the maximum absolute value of the response allows 


the excitation waveform to be positioned with the first point 
of the wire. 
2. Thin Wire Measured Data 

Measurements were performed in the Transient 
Electromagnetic Scattering Lab (TESL) to test the algorithm by 
using a thin wire with the same dimensions as the one used for 
the previous simulation. A detailed explanation of the 
procedure and the techniques used for measuring the scattering 


response from the thin wire, as well as other scale model 
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targets may be found in Morgan and McDaniel [Ref. 16] 
and Bresani [Ref. 17]. 

One measurement was available for each of the aspects 
at 30, 45 and 90 degrees. The scattering response waveforms 
obtained from the measurements are illustrated in Figures 24, 
26 and 28. These waveforms were compared with the calculated 
waveforms obtained from the TDIE program. It can be seen from 
the figures that the natural modes between the calculated and 
the measured waveforms do not have exactly the same 
frequencies. 

The spectrum of the measured response waveforms were 
also obtained, giving a distribution of energy within the 
bandwidth 1-12 GHz on eight frequencies, as shown in Table 


i Bed a 6 


Table III DISTRIBUTION OF ENERGY WITHIN THE SPECTRUM OF 
THIN WIRE 


Aspect Frequencies Frequencies Frequencies 
[degrees] [GHz] with [GHz] with [GHz] with 
most energy less energy no energy 


30 4.4 550" 1.6,8.6,10.2 11.6 
2 


3, 
ae 


45 1.6,5.6 7.2,8.6mll0c2> 
11.6 


wee: 3,5.6,8.6, 
10.2,11.6 
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Figures 25, 27 and 29 illustrate the spectrum of the measured 
thin wire scattering response for each aspect. 

When the measured thin wire response was processed, 
the points used to run the algorithm were from 180 to 200, 
starting at the initial excitation point, while the number of 
asking poles ranged from 30 to 60. The best results were 
obtained when the number of asking poles approached 40. The 
feed-forward order of the system was either the same as the 
number of asking poles, or was calculated by determining the 
early-time length. The scaling of the driving waveform and 
the positioning of the excitation with respect to the response 
was computed, as previously described. Figure 30 illustrates 
the extraction of the poles from the TDIE response, while 
Figure 31 illustrates the extraction of the poles from the 
measured waveform for the combined aspects. Figures 32-37 
illustrate the extraction of the poles from the measured 
waveform for each aspect separately. The pole results 
obtained using both the bias compensation method developed in 
this thesis and Cadzow-Solomon's bias compensation method have 
been plotted in Figures 32-37 along with the poles obtained 
without any bias compensation method, so that the results from 


the three types of methods may be compared. 
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Figure 26 Thin wire TDIE & Measured Response, 45 deg aspect 
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Figure 27 Spectrum of thin wire measured response, 45 deg 
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Figure 28 Thin wire TDIE & Measured Response, 90 deg aspect 
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Figure 33 
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Figure 34 
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Figure 35 Poles at z-plane, 30 deg aspect 
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Figure 36 Poles at z-plane, 45 deg aspect 


66 


+:TDIE response 


O:measured response 
no bias compensation 


x measured response 
Lazarakos’s bias compensation 


*:measured response 
Cadzow-Solomon’s bias compensation 





Figure 37 Poles at z-plane, 90 deg aspect 
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IV. POLES FROM SCALE MODEL AIRCRAFT 


Scattering data for four different scale model aircraft 
targets were processed without bias compensation, using the 
Cadzow-Solomon algorithn. Poles were extracted from the 
measured scattering responses of the aircraft targets to 
double-Gaussian electromagnetic excitation incident at 0, 30, 
90, and 180 degrees. The O degree aspect represents a nose-on 
measurement, while the 90 degree represents a broadside 
measurement. Signatures are shown in Figures 38-45 for targets 
1 and 2. 

The bandwidth of the TESL, 1-12 GHz, was matched to the 
scaling factor of the model aircraft targets. This scaling 
factor was 1/72 for all of the model aircraft targets used. 
Table IV contains the full scale, significant dimensions of 
these aircraft targets. The aircraft data was collected by 
Bresani (Ref. 17] at the four incident angles, as previously 


described. 


Table IV SCALED DIMENSIONS OF AIRCRAFT TARGETS 


Target number 


Overall length (cm) 
Wingspan (cm) 
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A. DATA PROCESSING 

The model aircraft scattering data was processed using a 
number of time points ranging from 200 to 240. The algorithm 
processing was started at the point of initial excitation. 

The positioning of the excitation waveform was set up 
manually with respect to the response waveform from observa- 
tions of the scattering data. Manual positioning for the 
driving waveform was required as no obvious condition existed 
for the model aircraft, as was the case for the thin-wire. 

The number of asking poles was set to 60, as the value of 
60 was found from previous experimentation to represent the 
most suitable results. The number of asking zeros was the 
same as the number of asking poles. This value was based on 
the fact that the computed early-time for the aircraft target 
is usually about 2L/c. Conversion of this early-time to the 
appropriate number of time points (equal to the number of 
asking zeros) provides a number larger than the number of 
asking poles. However, the number of asking zeros may not be 
larger than the number of asking poles, the two values are set 


to be equal. 


B. RESULTS FROM EXTRACTING THE POLES 
The model aircraft scattering data showed that the poles 
were less likely to group than the thin-wire data. The highly 


complex nature of the aircraft scatterer combined with the 


69 


small number of significant data points (about 150) became a 
difficult problem for the algorithm to solve. 

Figures 46-47 illustrate the spectrum of the data record 
for aircraft targets 1 and 2, for different incident angles. 

Figures 48-51 illustrate the extraction of the poles from 
the same aircraft target for the combined aspects. To obtain 
an initial indication of the possibility of target identifi- 
cation by pole extraction, nose-on measurement results have 
been compared in Figures 52-53. 

The thin wire data (Chapter III) indicated that the 
principle poles provided well defined clusters for all the 
incident angles examined, despite the use of a wide range of 
processing parameters. However, the scale model aircraft 
targets did not provide any well defined clusters under the 
limited attempts at pole estimation conducted here. The prin- 
Ciple reason for this difference appears to be that the scale 
models have more complicated pole patterns than the thin wire 
targets and insufficient time was available to explore ideas 
for optimizing the performance of the estimation algorithm. 

Processing the experimental scattering data for aircraft 
models, as performed herein, is only an initial attempt to 
demonstrate resonance estimation for real-world targets 
embodying complex configurations. Processing for this 
aircraft data was conducted for only one Wee at the end of 
the thesis research. It is recommended that this same data be 


used for a follow-on thesis. 
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Figure 51 Poles extracted from model aircraft 2 measured 
response in z-plane 
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V. SUMMARY AND CONCLUSIONS 


A. SUMMARY 

This thesis has attempted to demonstrate radar target 
identification by building on the earlier work of Norton [Ref. 
11), Larison (Ref. 12] and Murphy [Ref. 13]. The largest 
portion of the work consisted of testing the Cadzow-Solomon 
algorithm using synthetic and thin wire measurement data. 

Chapter I introduced the use of the Kumaresan-Tufts and 
the Cadzow-Solomon algorithms to locate target poles for a 
Non-Cooperative Target Recognition system. The resonance- 
based radar target identification problem was discussed and 
the transient electromagnetic scattering described. 

Chapter II consisted of two parts. The first part 
discussed early methods used to solve target identification 
problems. This discussion included a description of Prony's 
method and also, Singular Value Decomposition on which the 
Kumaresan-Tufts and Cadzow-Solomon algorithms were developed. 
The Kumaresan-Tufts algorithm was Wareiiaa ad in detail, 
including Norton's bias compensation method. The performance 
of the Kumaresan-Tufts algorithm was also demonstrated, as 
illustrated in Figures 3 through 10. The second part 


described the Cadzow-Solomon algorithm. Two bias compensation 
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methods were included: that of Cadzow-Solomon and that of the 
author. 

Chapter III demonstrated the Cadzow-Solomon algorithm 
testing in two phases. The first phase of testing was 
performed with synthetic data, while the second phase was 
performed with thin wire measurement data. The thin wire data 
testing phase attempted to evaluate conditions appearing from 
a real target response. Results of the synthetic signal model 
testing are illustrated in Figures 13 through 23. The results 
of the thin wire scattering signature testing are illustrated 
in Figures 24 through 37. 

Chapter IV considered an initial attempt at extraction of 
poles from scale model aircraft measurements obtained in the 
NPS Transient Electromagnetic Scattering Lab. Measurement 
data was processed using the Cadzow-Solomon algorithm and the 
unsuccessful results are illustrated in Figures 38 through 53. 

Testing of both the Kumaresan-Tufts and the Cadzow-Solomon 
algorithms was performed using MATLAB programs. A sequence of 
programs was written to complete the demonstration of the 
target identification problem. Appendices A to E present only 
a part of this sequence of programs, including the theoretical 


models of the Cadzow-Solomon algorithn. 
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B. CONCLUSIONS 

Both the Kumaresan-Tufts and the Cadzow-Solomon algorithms 
can effectively extract poles from the scattering response of 
Simple high-Q targets such as thin wires. As the late-time 
portion of the response signal is weak (with low SNR) and the 
Cadzow-Solomon algorithm has the ability to use the early-time 
portion, this thesis concentrated mainly on the use of the 
Cadzow-Solomon algorithn. 

The Cadzow-Solomon algorithm extracts poles of 
synthetically generated data, integral equation computed data, 
and thin wire measured scattering data with excellent results. 
A signal-to-noise ratio above 10 dB is required, depending on 
the damping rate of the data. The system order is 
intentionally overestimated. The excess order provides the 
flexibility to model the noise and improves the estimation of 
parameters of exponentially damped sinusoidal signals in 
noise. The Singular Value Decomposition method alleviates 
severe ill-conditioning of the data matrix. Backward 
prediction and SVD are used to separate the computed signal 
poles from spuriously computed noise poles introduced by the 
overestimated system order. 

The most critical parameters required for the successful 
thin-wire processing were the selection of the appropriate 
starting point to begin the data processing and the 
appropriate system order. The best results were obtained with 


the starting point of the data processing set at the 
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excitation starting point, and with the processing system 
order at about 3 times the true system order. 

The existence of noise in the target's response produces 
a bias in the positioning of the extracted poles. Thus, 
several bias compensation schemes may be developed. 

Given a very limited initial effort, the Cadzow-Solomon 
algorithm was unable to extract poles of scale model aircraft 
scattering data with adequate accuracy. It is conjectured 
that the data points in the natural mode information response 
are too few for the algorithm to model both the target poles 


and the noise poles. 
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APPENDIX A. ARMA COEFFICIENT AND RECURSIVE SIGNAL GENERATOR 


PROGRAM LISTING 
ARMAGEN1 .M 
Program generates a(k), b(k) coefficients of 
BO) +b (1) 427-14 2.24 (q) *2°-q 
ira (a) s2-—it. tat Pp) *Z°—p 
given the s-plane poles, residues, number of time points 


and time window. 
Programed by Gregory Lazarakos, 5 Mar 1991. 


GP JP IP dP OP dP OP dP OP dP OP dP 
e 
wm 
N 
ll 
' 
] 


'del temp.mat 

case='synt'; 
dt=t0./(notp-1); 
i=sqrt(-1); 
dm=exp ( (sigma+i*omega) *dt) ; 
alpha=ampl.*exp(i*phase) ; 
dm=dm'; 
alpha=alpha'; 
['Please wait...'] 
[b,a]=residue(alpha,dm,0) ; 
=real(a); 
B=real(b); 


% Program generates time response of an ARMA system 
% via the equation 
% N L 
* y(n) = SUM A(k) *y(n-k+1) + SUM B(k) *x(n-k+1) ,for n=1l:notp 
$ k=2 k=1 
% given A(k), B(k) and input excitation record x(n). 
$ Programed by Gregory Lazarakos, 10 Mar 1991. 

N=size(A) ; 

N=N (2) ? 

L=size(B) ; 

L=L(2); 

for j=2:N 

A(})=-A(j);? 

end 

% 


['Excitation is double Gaussian'] 
excgen 
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[{'Please wait...'] 


ye=zeros(1,notp); 
ye (1) =B(1) *x(1) ; 
energy=(ye(1)).*2;7 
for n=2:notp 
ye (n)=0.0; 
Ln=(n;Lj; 
Lmax=min (Ln) ; 
for k=1:Lmax 
ye (n)=ye(n)+B(k) *x(n-k+1) ; 
end 
Nn=({n;N]; 
Nmax=min(Nn) ; 
for k=2:Nmax 
ye (n) =ye (n) +A (k) *ye (n-k+1) ; 
end 
energy=energy+ (ye(n)).%*2; 
end 
rms=sqrt(energy) ; 
i1c=zeros(1,notp); 
for n=l1:notp 
ic(n)=ye(n)./rms; 


en 
% 
axis([O notp -0.4 0.8]) 
plot(ic) 
title([{'Generated signal w/o noise, medium Q, by 


',int2str(nop),' poles']); 
Xlabel('time points'); 
ylabel('amplitude rms'); 


op ap Re i 
pause 
hard=input('Do you want a hardcopy for the plot? [n] : 
ie eSote iee 
if hard=='y' 
plotfile=input('Enter filename : ','s'); 
eval([({'meta ',plotfile]) 
end 
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APPENDIX B. DOUBLE GAUSSIAN SMOOTHING FUNCTION GENERATOR 


PROGRAM LISTING 


WP dP dO WO dP AO AO AP AO AP OP OP AO OP 


af 


EXCGEN.M 


Program generates excitation record x(n) for a response 
from TESL data and synthetic data, 
via the equation of a double Gaussian waveform 


x(n) = Al*exp(-al*t*2) - A2*exp(-a2*t*2) 
given the time window. 


Program uses input values of the 10% height of the low 
and high frequency ends of the double Gaussian frequency 


response to determine the pulse widths in the time domain. 


Programed by Gregory Lazarakos, 8 Apr 1991. 
Last Revision August 6, 1991. 
case=='meas' 
if filename(1)=='F' 
LENGTH=0.25; 
elseif filename(siz-1:siz)=='sp' 
LENGTH=eval (filename (3:4)); 
LENGTH=LENGTH/100; 
else 
LENGTH=0.1; 
end 


end 


npts=notp; 

tmin=0.0; 

tmax=to; 

DT=dt; 

T1=0.147; 

T2=0.314; 

thr=input('Enter threshold value for time in nsec 


[1.25]:'); 


gg 


if isempty (thr) 
thr=1.25; 
end 
al=(4.0.*10g(10))./(T1.%2); 
a2=(4.0.*10g(10))./(T2.%2); 
Al=sqrt(al)./(sqrt(al)-sqrt(a2)); 
A2=A1-1.0; 
case=='synt' 
point=input('Which to be the point, the excitation 


starts?'); 


X=zeros(1,npts) ; 
for n=l1:npts 
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t=(n-11-point) .*DT; 
if t<=thr 
€1=Al.*exp(-al.*t.%2) ; 
e2=A2.*exp(-a2.*t.%2) ; 
x (n) =el-e2; 
else 
x(n) =0.0; 
end 
end 
elseif case=='meas' 
asprad=aspect* (pi/180) ; 
top=input('Do you want to position the excitation, auto 
manual ¢f{a,m) ; 9°, Sie, 
if top=='a' 
[maxic, index1]=max(ic) ; 
[minic, index2]=min(ic) ; 
X=zeros(1,npts) ; 
tdel=cos (asprad) *(2*LENGTH/3e8) ; 
ndel=round (tdel/ (dt*le-9) ); 
if abs(maxic) >=abs(minic) 
point=index1-11-ndel; 
for n=1:npts 
t= (n-indexi+ndel) .*DT; 
if abs(t) <=10*DT 
e1l=Al.*exp(-al.*t.%*2); 
e2=A2.*exp(-a2.*t.%2) ; 
Xx (n) =el-e2; 
else 
x(n) =0.0; 
end 
end 
X=x*0.6; 
else 
point=index2-11-ndel; 
for n=l1:npts 
t=(n-index2+ndel) .*DT; 
if abs(t)<=10*DT 
el=Al.*exp(-al.*t.%2) ; 
e2=A2.*exp(-a2.*t.%*2) ; 
x (n) =el-e2; 
else 
x(n)=0.0; 
end 
end 
X=x*0.6; 
end 
end 
if top=='m!' 


or 


point=input('Which to be the point, the excitation 


starts 7 "%); 
X=zeros(1,npts) ; 
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for n=l:npts 
t=(n-11-point) .*DT; 
if abs(t)<=10*DT 
e1=Al.*exp(-al.*t.%2); 
e2=A2.*exp(-a2.*t.%2) ; 
x(n) =el-e2; 
else 
x(n)=0.0; 
end 
end 
X=X*.6; 
end 
ltst=(1+cos(aspraqd) ) *(LENGTH/3e8) ; 
nltst=round(ltst/ (dt*le-9) ); 
ltimp=point+22+nltst; 
else 
exit 
end 


if case=='synt'! 
axis([{[O notp -0.4 1}) 
plot(x) 
title('Double-Gaussian Smoothing Function'); 
xlabel('time points'); 
ylabel('amplitude rms'); 
elseif case=='meas' 
pLOtiemDes, x, '——",1:notp,ic,1ltimp,ic(ltimp),'*') 
if filename(1)=='F' 
title(['Data response signal from aircraft',filename(1:3),', 
aspect=!',num2str(aspect),' deg ')); 


elseif filename(siz-1:siz)=='sp' 
title(['Data response signal from sphere', filename(3:4),'cm, 
aspect=',num2str(aspect),' deg ']); 
else 
title([{'Data response Signal from thin wire, 
aspect=!',num2str(aspect),' deg ']); 
end 
xlabel ('solid=response dash=excitation')j; 
ylabel('amplitude '); 
text (ltimp, 0.1, '----- late time----- >"); 


if filename(1:2)=='et' 
text(100,0.4, 'TDIE response'); 
end 
grid; 
else 
exit 
end 
pause 
hard=input('Do you want a hardcopy for the plot? 
[n]:','s'); 
if hard=='y' 
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plotfile=input('Enter filename : ','s'); 
eval(['meta ',plotfile}) 
end 
if case=='meas' 
Save temp ic energy x point dt tO notp filename aspect siz 
LENGTH 
end 
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APPENDIX C. SYNTHETIC NOISE GENERATOR 


PROGRAM LISTING 
NOISEGEN .M 


Synthetic Noise Generation 
Programmed by Greg Lazarakos 20/2/91 


WH AP AP HO AO 


SN=input('How many db ? (7,10,20,30) [20]:'); 
if isempty (SN) 
SN=20; 
end 
sn=SN./10; 
sn=10.%sn; 
en=0; 
for k=l:notp 
en=en+(ic(k)).*2; 
end 
ava=en./(notp.*sn) ; 
ava=sqrt (ava) ; 
rand('normal') 
w=rand(ic) *ava; 
ic=ic+w; 
% 
save temp ic x point filename w ye energy dt dm notp omega 
sigma to 
3 : 
axis([O notp -0.4 0.8}) 
plot(ic) 
title({'Generated signal with noise S/N=',int2str(SN),' db']}); 
xlabel('time points'); 
ylabel('amplitude rms'); 
grid; 
pause 
hard=input('Do you want a hardcopy for the plot? 
iis ', 's'); 
if hard=='y' . 
plotfile=input('Enter filename : ','s'); 
eval(['meta ',plotfile}) . 
end 
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APPENDIX D. SPECTRUM ESTIMATION 


PROGRAM LISTING 


% 
% 
% 
% 


SPECTRUM .M 


Programed by Gregory Lazarakos, 16 Apr 1991. 
Last revision, 6 Aug 1991. 
df=1./((notp-1) .*dt) ; 

Ic=fft(ic,notp) ; 

f=(O:notp-1) *df; 

fmax=1./dt; 

axis([{[0O fmax -3 3)) 

plot (f£ 1¢) 


title('Fourier transform of the response signal'); 


xlabel('Frequency in GHz'); 
ylabel('Amplitude'); 

grid; 

pause 

SIC=(abs(IC)).*2; 

SIC=SIC/notp; 

axis({[0 20 0 max(SIC) }) 

DLOoe(L sie) 

if isempty (SN) 

title('Spectrum of the response signal '); 
else 

title(['Spectrum OF the response signal 


S/N=',int2str(SN),' db'}); 


(n] 


end 
Xlabel('Frequency in GHz'); 
ylabel('Amplitude') ; 
grid; 
pause 
hard=input('Do you want a hardcopy for the plot? 
. re ae 
if hard=='y' 
plotfile=input('Enter filename : ','s'); 
eval([{'meta ',plotfile}) 
end 
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with 


APPENDIX E. THE CADZOW-SOLOMON POLE EXTRACTION ALGORITHM 


PROGRAM LISTING 


JP dP dP dP dO AP dP OP 


CADSOL1 .M 


Cadzow-Solomon's Algorithm for Extracting the Poles 
and Residues 

Version 1.0 

Programmed by Greg Lazarakos 3/16/91 


Set the first sample (kapa) 
if point>0 

kapa=point; 
else 

kapa=1; 
end 
Set the number of samples (CN) 
CN=input('How many samples to process ? (>120) [200]:'); 
1f isempty(CN) 

CN=200; 
end 
Compute the SNR for the processed time window 
en=0; 
nen=0; 
for k=kapa:CN 

en=en+(ic(k)).%*2; 

nen=nen+ (w(K) ).*2; 
end 
SNR=en./nen; 
SNR=10.*10g10(SNR) ; 
Set the number of poles (CL) > number of real poles (CM) 
CM <= CL and 2*CL <= CN-CL 
CL=input('How many unknown poles ? (>20) [60}: '); 
if isempty(CL) 

CL=60; 
end 
Set the number of expected zeros (CK) : 
CK=input('How many expected zeros ? [No. poles]: '); 
if isempty (CK) 

CK=CL 
end 
('Please wait...'] 
Forming the matrix YI[(CN-CL) *CL] from the data ic(n) 

YI=zeros(CN-CL,CL) ; 
for 4=kapa+1:CN-CLtkapa 

YI(j-kKapa,:)=ic(j:j+CL=-1) ; 
end 
Forming the matrix h[(CN-CL)*1]} from the data ic(n) 
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% 
% 
% 


h=[( J; 

h=ic(kapa:CN-CLt+kapa-1) ; 

h=h'; 

Forming the matrix XI[(CN-CL)*CK] from the data x(n) 
XI=zeros (CN-CL, CK) ; 

for j=1:CN-CL 
XI(j,:2)=x(j:j+CK-1) ; 

end 

Unify matrices YI and XI as I=[YI!xI]} 

I=[YI XI}; 

Set the tolerance 

tol=0.000001; 

Find the vector of backward prediction coefficients 

beta=pinv(I,tol)*([h]; 

Set the coefficients of the prediction-error filter 


polynomial 


ca=zeros(1,CLt1l); 
ca(1)=1; 
for j=1:CL 
ca(jt+1)=-beta(j); 
end 
Set the coefficients of the polynomial B(z) 
cb=zeros(1,CK); 
for j=1:CK 
cb(j})=beta(j+CL) ; 
end 
Find the residues and the poles 
[resid,d,gk]=residue(cb,ca); 
Find the signal poles 
S=log(d)/dt; 


1=0; 
for j=1:CL 
if real(s(j))>0, 
k=k+1; 
polesi1(k)=-s(j)? 
residi(k)=resid(j); 
else 
1=1+1; 
nopoli1(1)=-s(3);3 
noresi(1)=resid(j) ; 
end 
end 


dmi=exp (polesi*dt) ; 

dm2=exp (nopoli*dt) ; 

alphal=exp(residi*dt) ; 

CM=k ; 

['The signal has the following ',num2str(CM),' real poles.'] 
polesil 
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